// Copyright (c) Lawrence Livermore National Security, LLC and other VisIt
// Project developers.  See the top-level LICENSE file for dates and other
// details.  No copyright assignment is required to contribute to VisIt.

// ************************************************************************* //
//                             avtEulerianQuery.C                           //
// ************************************************************************* //

#include <avtEulerianQuery.h>

#include <set>
#include <vector>

#include <vtkCell.h>
#include <vtkDataSet.h>
#include <vtkDataSetRemoveGhostCells.h>
#include <vtkGeometryFilter.h>
#include <vtkPolyData.h>

#include <avtQueryableSource.h>

#include <DebugStream.h>
#include <ImproperUseException.h>
#include <NonQueryableInputException.h>

#ifdef PARALLEL
#include <mpi.h>
#include <avtParallel.h>
#endif

using std::set;
using std::vector;
using std::string;



// ****************************************************************************
//  Method: avtEulerianQuery::avtEulerianQuery
//
//  Purpose:
//      Construct an avtEulerianQuery object.
//
//  Programmer:   Akira Haddox
//  Creation:     June 28, 2002
//
//  Modifications:
//    Kathleen Bonnell, Fri Nov 15 12:37:11 PST 2002 
//    Moved from avtEulerianFilter.  
//
// ****************************************************************************

avtEulerianQuery::avtEulerianQuery()
{
    gFilter = vtkGeometryFilter::New();
}

// ****************************************************************************
//  Method: avtEulerianQuery::~avtEulerianQuery
//
//  Purpose:
//      Destruct an avtEulerianQuery object.
//
//  Programmer:   Akira Haddox
//  Creation:     June 28, 2002
//
//  Modifications:
//    Kathleen Bonnell, Fri Nov 15 12:37:11 PST 2002 
//    Moved from avtEulerianFilter. 
//
// ****************************************************************************

avtEulerianQuery::~avtEulerianQuery()
{
    if (gFilter)
        gFilter->Delete();
}


// ****************************************************************************
//  Method: avtEulerianQuery::VerifyInput
//
//  Purpose:
//    Rejects non-queryable input && input that has topological dimension == 0
//
//  Programmer: Kathleen Bonnell
//  Creation:   September 3, 2004
//
//  Modifications:
//
// ****************************************************************************

void
avtEulerianQuery::VerifyInput()
{
    avtDataObjectQuery::VerifyInput();
    if (GetInput()->GetInfo().GetAttributes().GetTopologicalDimension() == 0)
    {
        EXCEPTION1(NonQueryableInputException,
            "Requires plot with topological dimension > 0.");
    }
}





// ****************************************************************************
//  Method: avtEulerianQuery::Execute
//
//  Purpose:
//      Computes the Eulerian number of the input dataset. 
//
//  Arguments:
//      inDS      The input dataset.
//      dom       The domain number.
//
//  Programmer:   Akira Haddox
//  Creation:     June 28, 2002
//
//  Modifications:
//    Kathleen Bonnell, Fri Nov 15 12:37:11 PST 2002 
//    Moved from avtEulerianFilter::DeriveVariable.
//    Instead of creating a data array, keep track of the Eulerian
//    for this domain in domToEulerMap. 
//
//    Hank Childs, Tue Aug 24 08:47:57 PDT 2004
//    Manually remove ghost cells, since the geometry filter only removes
//    vtkGhostLevels, not vtkGhostNodes.
//
//    Kathleen Bonnell, Mon May 16 10:28:15 PDT 2005 
//    Fix memory leak. 
//
//    Kathleen Bonnell, Wed May 17 15:22:06 PDT 2006 
//    Remove call to SetSource(NULL) as it now removes information necessary
//    for the dataset. 
//
//    Kathleen Biagas, Fri Jan 25 16:35:54 PST 2013
//    Call Update on filter, not data object.
//
// ****************************************************************************

void 
avtEulerianQuery::Execute(vtkDataSet *in_ds, const int dom)
{
    vtkDataSet *nds = (vtkDataSet*) in_ds->NewInstance();
    nds->ShallowCopy(in_ds);
    //nds->SetSource(NULL);

    gFilter->SetInputData(nds);
    vtkDataSetRemoveGhostCells *ghost_remover =
                                             vtkDataSetRemoveGhostCells::New();
    ghost_remover->SetInputConnection(gFilter->GetOutputPort());
    ghost_remover->Update();
    vtkDataSet *out = ghost_remover->GetOutput();
    if (out->GetDataObjectType() != VTK_POLY_DATA)
    {
        // The output of a geometry filter should always be poly-data.
        EXCEPTION0(ImproperUseException);
    }

    vtkPolyData *pds = (vtkPolyData *) out;

    // I believe this isn't good enough. I believe the facelist filter
    // simply passes points through, and only modifies the cell structure.
    // As such, the number of points is greater than the number of points
    // that are actually used.
    int nPoints = pds->GetNumberOfPoints();
    // So we'll keep track manually of which points are used.
    vector<bool> pointsUsed(nPoints,false);

    int nCells = pds->GetNumberOfCells();

    int Eulerian;

    // We go through all the cells, and for each cell, add each edge of
    // that cell into our set of edges. Due to the nature of a set, and
    // that the ordering of nodes in an edgepair not mattering was specified
    // in the edgepair class, we will basically be able to tell how many
    // edges are in the mesh simply by the size of the set.

    set<edgepair> edges;

    for (size_t i = 0; i < (size_t)nCells; i++)
    {
        vtkCell *cell=pds->GetCell(i);
        int numCellPoints = cell->GetNumberOfPoints();

        // Marked all the points in the cell as being used in our records
        for (int m = 0; m < numCellPoints; m++)
            pointsUsed[cell->GetPointId(m)] = true;

        switch (cell->GetCellType())
        {
        // Triangles, straight forward. Three points, three edges
            case VTK_TRIANGLE:
                edges.insert(edgepair(cell->GetPointId(0),
                                      cell->GetPointId(1)));
                edges.insert(edgepair(cell->GetPointId(1),
                                      cell->GetPointId(2)));
                edges.insert(edgepair(cell->GetPointId(2),
                                      cell->GetPointId(0)));
                break;

        // Quads, four points, four edges.
            case VTK_QUAD:
                edges.insert(edgepair(cell->GetPointId(0),
                                      cell->GetPointId(1)));
                edges.insert(edgepair(cell->GetPointId(1),
                                      cell->GetPointId(2)));
                edges.insert(edgepair(cell->GetPointId(2),
                                      cell->GetPointId(3)));
                edges.insert(edgepair(cell->GetPointId(3),
                                      cell->GetPointId(0)));
                break;

        // Polygons are also easy. There's an edge between node 0 and N,
        // and all the other edges are bewtween edge n and n+1
            case VTK_POLYGON:
                edges.insert(edgepair(cell->GetPointId(0),
                                      cell->GetPointId(numCellPoints-1)));
                for (size_t j = 1; j < (size_t)numCellPoints-1; j++)
                    edges.insert(edgepair(cell->GetPointId(j),
                                          cell->GetPointId(j+1)));
                break;

        //Each node [save for outer ones] is connected to two nodes
        //before them, and two after. Since we don't repeat edges,
        //we can just count edges infront of us to get them all
            case VTK_TRIANGLE_STRIP:
                for (size_t j = 0; j < (size_t)numCellPoints-2; j++)
                {
                    edges.insert(edgepair(cell->GetPointId(j),
                                          cell->GetPointId(j+1)));
                    edges.insert(edgepair(cell->GetPointId(j),
                                          cell->GetPointId(j+2)));
                }
                break;

            default:
                debug5 << "Geometry filter returned unexpected type." << endl;
                debug5 << "Type is: " << cell->GetCellType() << endl;
        }
    }    

    int numUsedPoints = 0;
    for (size_t i = 0; i < pointsUsed.size(); i++)
        if (pointsUsed[i])
            ++numUsedPoints;

    // Now for the magic: Euler-Descartes formula
    Eulerian = numUsedPoints - (int)edges.size() + nCells;

    domToEulerMap.insert(DomainToEulerMap::value_type(dom, Eulerian));
    ghost_remover->Delete();
    nds->Delete();
}



// ****************************************************************************
//  Method: avtEulerianQuery::PostExecute
//
//  Purpose:
//      This is called after all of the domains are executed.
//      If in parallel, gathers the Eulerian number for each domain onto one 
//      processor. 
//      Constructs an output string listing Eulerian number for each domain, 
//      for output to the gui window.
//
//  Programmer: Kathleen Bonnell
//  Creation:   November 15, 2002
//
//  Modifications:
//    Kathleen Bonnell, Fri Jul 11 16:31:45 PDT 2003
//    Renamed 'SetMessage' to 'SetResultMessage'.
//
//    Kathleen Bonnell, Mon Dec 22 16:45:56 PST 2003 
//    Use retrieved DomainName if available. 
//
//    Mark C. Miller, Wed Jun  9 21:50:12 PDT 2004
//    Eliminated use of MPI_ANY_TAG and modified to use GetUniqueMessageTags
//
//    Mark C. Miller, Mon Jan 22 22:09:01 PST 2007
//    Changed MPI_COMM_WORLD to VISIT_MPI_COMM
//
//    Kathleen Biagas, Thu Feb 13 15:04:58 PST 2014
//    Add Xml results.
//
// ****************************************************************************

void
avtEulerianQuery::PostExecute(void)
{
#ifdef PARALLEL
    int myRank, numProcs;
    int size, i, j;
    int *buf;
 
    MPI_Comm_rank(VISIT_MPI_COMM, &myRank);
    MPI_Comm_size(VISIT_MPI_COMM, &numProcs);
    int tags[2];
    GetUniqueMessageTags(tags, 2);
    int mpiSizeTag = tags[0];
    int mpiDataTag = tags[1];
    if (myRank == 0)
    {
        for (i = 1; i < numProcs; i++)
        {
            MPI_Status stat, stat2;
             
            MPI_Recv(&size, 1, MPI_INT, MPI_ANY_SOURCE,
                     mpiSizeTag, VISIT_MPI_COMM, &stat);
            if (size > 0)
            {
                buf = new int [size];
                MPI_Recv(buf, size, MPI_INT, stat.MPI_SOURCE, mpiDataTag,
                         VISIT_MPI_COMM, &stat2);
                for (j = 0; j < size/2; j++)
                { 
                    domToEulerMap.insert(DomainToEulerMap::value_type(
                        buf[j*2], buf[j*2+1]));
                } 
                delete [] buf;
            }
        }
    }
    else
    {
        size = domToEulerMap.size() * 2;
        MPI_Send(&size, 1, MPI_INT, 0, mpiSizeTag, VISIT_MPI_COMM);
        if (size > 0)
        {
            buf = new int[size];
            DomainToEulerMap::iterator iter;
            for (iter = domToEulerMap.begin(), i = 0; 
                 iter != domToEulerMap.end(); iter++, i++)
            {
                buf[2*i] = (*iter).first;
                buf[2*i+1] = (*iter).second;
            }
            MPI_Send(buf, size, MPI_INT, 0, mpiDataTag, VISIT_MPI_COMM);
            delete [] buf;
        }
        return;
    }
#endif

    string msg;
    char msgBuff[500];
    DomainToEulerMap::iterator iter;
    int blockOrigin = GetInput()->GetInfo().GetAttributes().GetBlockOrigin();
    MapNode de;
    for (iter = domToEulerMap.begin(); iter != domToEulerMap.end(); iter++)
    {
        string domainName;
        std::stringstream dns;
        GetInput()->GetQueryableSource()->GetDomainName(
            queryAtts.GetVariables()[0], queryAtts.GetTimeStep(), 
            (*iter).first, domainName);
        if (domainName.size() > 0)
        {
            snprintf(msgBuff, 500, "Eulerian for %s is %d\n", 
                     domainName.c_str(), (*iter).second);
            dns << domainName;
        }
        else
        {
            snprintf(msgBuff, 500, "Eulerian for domain %d is %d\n", 
                     (*iter).first + blockOrigin, (*iter).second);
            dns << "domain " << (*iter).first+blockOrigin;
        }
        de[dns.str()] = (*iter).second;
        msg += msgBuff;
    }
    if (msg.size() == 0)
        msg = "Eulerian could not compute.\n" ;
    SetResultMessage(msg);
    if (de.GetNumEntries() > 0)
    {
        MapNode result_node;
        result_node["eulerian"] = de;
        SetXmlResult(result_node.ToXML());
    }
}


// ****************************************************************************
//  Method: avtEulerianQuery::PreExecute
//
//  Purpose:
//      This is called before any of the domains are executed.
//      Ensures the storage for Eulerian numbers is clear. 
//
//  Programmer: Kathleen Bonnell
//  Creation:   November 15, 2002
//
//  Modifications:
//    Jeremy Meredith, Thu Feb 15 11:55:03 EST 2007
//    Call inherited PreExecute before everything else.
//
// ****************************************************************************

void
avtEulerianQuery::PreExecute()
{
    avtDatasetQuery::PreExecute();

    if (!domToEulerMap.empty())
        domToEulerMap.clear();
}


